Step 1: Problem Formulation¶
To formulate the problem, we first have to import the dike model. This file is named 'dike_model_function' and this simulates the whole chain of functions and files. This basically is the realtionships in the system (R).
Secondly, the file 'problem_formulation_group6' is imported, in which the problem is formulated. It prepares the uncertainties (X), levers (L) and performance metrics (M). The table below specifies the ranges for the deeply uncertain factors:
| Uncertainties | Range | Unit |
|---|---|---|
| Maximum breach height per location (Bmax) | 30 – 350 | m |
| Breach growth rate per location (Brate) | (1, 1,5, 10) | m |
| Discount rate | (1.5, 2.5, 3.5, 4.5) | |
| Probability of dike failure per location (pfail) | 0 – 1 | |
| Flood wave shape | 0 – 140 |
The table below specifies the levers of the model:
| Levers | Range |
|---|---|
| Room for the River per location (RfR) | 0-1 |
| Dike heightening per location | 0-10 |
| Early warning system (EWS) | 0-4 |
Within this file, the performance metrics is adjusted to the prefered outcomes of interests of this report. A 7-objective problem formulation is created, which gives different costs as output as well as the expected number of deaths for the geographical areas: Upstream Gelderland, Downstream Gelderland and Overijssel. The table below shows all objectives within the problem formulation:
| Performance metrics |
|---|
| Expected annual damage |
| Dike investment costs |
| RfR investment costs |
| Evacuation costs |
| Expected annual deaths per geographical area |
Importing relevant packages and initialise the model¶
In [1]:
# Import the relevant packages
import numpy as np
import scipy as sp
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import networkx as nx
# Set a coloscheme to generate consistent looking plots
colorscheme = "viridis"
In [2]:
# Import the EMA workbench packages needed to specify the problem formulation in this notebook. Also import the dike_model_function
from ema_workbench import (Model, CategoricalParameter,
ScalarOutcome, IntegerParameter, RealParameter)
from dike_model_function import DikeNetwork # @UnresolvedImport
def sum_over(*args):
return sum(args)
In [3]:
# Import other relevant EMA workbench functions
from ema_workbench import (Model, MultiprocessingEvaluator, Policy, Scenario)
from ema_workbench.em_framework.evaluators import perform_experiments
from ema_workbench.em_framework.samplers import sample_uncertainties
from ema_workbench.util import ema_logging
import time
from problem_formulation_group6 import get_model_for_problem_formulation
ema_logging.log_to_stderr(ema_logging.INFO)
#choose problem formulation number, between 0-5
#each problem formulation has its own list of outcomes
dike_model, planning_steps = get_model_for_problem_formulation(2)
In [4]:
#enlisting uncertainties, their types (RealParameter/IntegerParameter/CategoricalParameter), lower boundary, and upper boundary
for unc in dike_model.uncertainties:
print(repr(unc))
uncertainties = dike_model.uncertainties
import copy
uncertainties = copy.deepcopy(dike_model.uncertainties)
In [5]:
#enlisting policy levers, their types (RealParameter/IntegerParameter), lower boundary, and upper boundary
for policy in dike_model.levers:
print(repr(policy))
levers = dike_model.levers
import copy
levers = copy.deepcopy(dike_model.levers)
In [6]:
#enlisting outcomes
for outcome in dike_model.outcomes:
print(repr(outcome))
Running the model through EMA workbench¶
Running the model through the EMA workbench using SequentialEvaluator. This code creates just 2 experiments, but that is fine. The results are not used for analysis.
In [7]:
#running the model through EMA workbench
from ema_workbench import (MultiprocessingEvaluator, ema_logging,
perform_experiments, SequentialEvaluator)
ema_logging.log_to_stderr(ema_logging.INFO)
with SequentialEvaluator(dike_model) as evaluator:
results = evaluator.perform_experiments(scenarios=1, policies=2)
In [8]:
#observing the simulation runs
experiments, outcomes = results
print(outcomes.keys())
experiments
Out[8]:
Running and exploring Base case scenario¶
Now we are going to run the base case scenario, in which no policy measures are implemented. This is run for 5000 scenarios.
In [10]:
from ema_workbench import Policy
policies = [Policy('No_measures', **{'0_RfR 0':0, '0_RfR 1':0, '0_RfR 2':0, 'A.1_DikeIncrease 0':0, 'A.1_DikeIncrease 1':0, 'A.1_DikeIncrease 2':0,
'1_RfR 0':0, '1_RfR 1':0, '1_RfR 2':0, 'A.2_DikeIncrease 0':0, 'A.2_DikeIncrease 1':0, 'A.2_DikeIncrease 2':0,
'2_RfR 0':0, '2_RfR 1':0, '2_RfR 2':0, 'A.3_DikeIncrease 0':0, 'A.3_DikeIncrease 1':0, 'A.3_DikeIncrease 2':0,
'3_RfR 0':0, '3_RfR 1':0, '3_RfR 2':0, 'A.4_DikeIncrease 0':0, 'A.4_DikeIncrease 1':0, 'A.4_DikeIncrease 2':0,
'4_RfR 0':0, '4_RfR 1':0, '4_RfR 2':0, 'A.5_DikeIncrease 0':0, 'A.5_DikeIncrease 1':0, 'A.5_DikeIncrease 2':0, 'EWS_DaysToThreat':0})]
In [11]:
#pass the policies list to EMA workbench experiment runs
from ema_workbench import (MultiprocessingEvaluator, ema_logging,
perform_experiments, SequentialEvaluator)
ema_logging.log_to_stderr(ema_logging.INFO)
n_scenarios = 5000
with MultiprocessingEvaluator(dike_model) as evaluator:
results = evaluator.perform_experiments(n_scenarios,
policies)
In [24]:
#observing the simulation runs
experiments, outcomes = results
experiments
# experiments
# Statistical analysis for all policies
# pd.DataFrame(outcomes).describe()
# Statistical analysis per policy
outcomes_df = pd.DataFrame.from_dict(outcomes)
outcomes_df = pd.concat([experiments['policy'], outcomes_df], axis=1)
df_step1 = pd.DataFrame(outcomes_df).describe()
In [27]:
df_step1
Out[27]:
The results of running the base case over 5000 scenario's shows that the number of casualties for three KPI's concerning expected deaths is way above the threshold set by the government. This implies that inspecting policy implementations is a viable thing to do.
In [25]:
df_step1.to_csv('step1_basecase_outcomes.csv')
In [26]:
total_df = pd.concat([experiments, outcomes_df], axis=1)
total_df
Out[26]:
In [28]:
total_df.to_csv('step1_basecase_outcomes_experiments.csv')
In [31]:
data = pd.DataFrame(outcomes)
fig, axs = plt.subplots(ncols=4, figsize = (20,6)) # Add sharey=True to have the same axis, makes more easy to visualize differences
sns.violinplot(data=outcomes_df, y='Expected Annual Damage', ax=axs[0])
sns.violinplot(data=outcomes_df, y='Expected Number of Deaths Upstream Gelderland', ax=axs[1])
sns.violinplot(data=outcomes_df, y='Expected Number of Deaths Downstream Gelderland', ax=axs[2])
sns.violinplot(data=outcomes_df, y='Expected Number of Deaths Overijssel', ax=axs[3])
plt.savefig('basecase.png')
plt.show()
As a more visual approach to analyzing the outcomes of interests, these violinplots were generated. These plots clearly indicate the distribution of the scenarios, in which is directly visible that the base case scenario results have a good possibility of a highly undesirable future situation.
To conclude, analyzing the base case scenario, in which no policy measures are implemented, provides clarity about the seriousness of the situation and the need for better flood risk management. It also provides insight into the current safety of the various areas and this, in combination with the improvements of the found policy implementations, can be useful at a later stage in a fair distribution of costs.
In [ ]:
# Saving the results of Step 1 for the base case scenario
outcomes_df.to_csv('step1_basecase_outcomes.csv')
Step 2: Searching for candidate solutions¶
In the second step of MORDM, candidate strategies are identified which are pareto optimal conditional on a reference scenario. These candiate strategies are identified through search with multi-objective evolutionary algorithms, that iteratively evaluate a large number of alternatives on multiple objectives until they find the best candidates.
In [12]:
from ema_workbench.em_framework.optimization import (EpsilonProgress)
convergence_metrics = [EpsilonProgress()]
In [13]:
start = time.time()
with MultiprocessingEvaluator(dike_model) as evaluator:
results2, convergence = evaluator.optimize(nfe=5e3, searchover='levers',
convergence=convergence_metrics,
epsilons=[1000,1000,1000,10,0.00001,0.00001,0.00001]
)
end = time.time()
print('The hypervolume optimization took ' + str(round((end - start)/60)) + ' minutes')
In [15]:
from ema_workbench.analysis import parcoords
outcomes = results2.loc[:, ['Expected Annual Damage', 'Dike Investment Costs', 'RfR Investment Costs', 'Evacuation Costs',
'Expected Number of Deaths Upstream Gelderland', 'Expected Number of Deaths Downstream Gelderland', 'Expected Number of Deaths Overijssel']]
limits = parcoords.get_limits(outcomes)
axes = parcoords.ParallelAxes(limits)
axes.plot(outcomes)
plt.show()
In [16]:
fig, (ax1) = plt.subplots(ncols=1, sharex=True, figsize=(8,4))
ax1.plot(convergence.nfe, convergence.epsilon_progress)
ax1.set_ylabel('$\epsilon$-progress')
ax1.set_xlabel('number of function evaluations')
plt.show()
In [18]:
# Save all pareto optimal solutions found
r = results2.iloc[:,7:]
r.to_csv('all_pareto_solutions.csv')
In [42]:
# Set constraints on the number of deaths and total costs to find optimal solutions
expected_number_of_deaths_constraint = 0.01
damage_costs_constraint = 5e8
d = results2.loc[:, 'Expected Number of Deaths Upstream Gelderland'].values + results2.loc[:, 'Expected Number of Deaths Downstream Gelderland'].values + results2.loc[:, 'Expected Number of Deaths Overijssel'].values
d = pd.DataFrame(d)
d = d[d[0] < expected_number_of_deaths_constraint]
r_1 = results2.loc[d.index.values.tolist(),:]
r_1['Total Investment Costs'] = r_1['Dike Investment Costs'] + r_1['RfR Investment Costs']
r_1 = r_1.drop(r_1[r_1['Total Investment Costs'] > damage_costs_constraint].index)
r_1
Out[42]:
In [43]:
r_1.to_csv('generated_policies.csv')
In [44]:
df = pd.read_csv('generated_policies.csv')
In [45]:
df2 = pd.DataFrame()
In [46]:
df2['Policy'] = df['Unnamed: 0']
df2['Expected Annual Damage'] = df['Expected Annual Damage']
df2['Dike Investment Costs'] = df['Dike Investment Costs']
df2['RfR Investment Costs'] = df['RfR Investment Costs']
df2['Evacuation Costs'] = df['Evacuation Costs']
df2['Expected Number of Deaths Upstream Gelderland'] = df['Expected Number of Deaths Upstream Gelderland']
df2['Expected Number of Deaths Downstream Gelderland'] = df['Expected Number of Deaths Downstream Gelderland']
df2['Expected Number of Deaths Overijssel'] = df['Expected Number of Deaths Overijssel']
In [48]:
outcomes = df2.loc[:, ['Expected Annual Damage', 'Dike Investment Costs', 'RfR Investment Costs', 'Evacuation Costs',
'Expected Number of Deaths Upstream Gelderland', 'Expected Number of Deaths Downstream Gelderland', 'Expected Number of Deaths Overijssel']]
limits = parcoords.get_limits(outcomes)
axes = parcoords.ParallelAxes(limits)
axes.plot(outcomes)
plt.show()
Step 3: Uncertainty Analysis¶
In [66]:
r_1 = pd.read_csv('generated_policies.csv')
In [67]:
r_2 = r_1.drop([o.name for o in dike_model.outcomes], axis=1)
r_2 = r_2.drop('Total Investment Costs', 1)
r_2 = r_2.rename(columns={'Unnamed: 0': 'Policy Number'})
r_2
Out[67]:
In [68]:
r_2.to_csv('policies_selection_step3.csv')
In [69]:
from ema_workbench import Policy
policies_to_evaluate = []
for i, policy in r_2.iterrows():
policies_to_evaluate.append(Policy(str(i), **policy.to_dict()))
In [70]:
n_scenarios = 1000
with MultiprocessingEvaluator(dike_model) as evaluator:
results3 = evaluator.perform_experiments(n_scenarios,
policies_to_evaluate)
In [71]:
experiments3, outcomes3 = results3
experiments3
# experiments
# Statistical analysis for all policies
# pd.DataFrame(outcomes).describe()
# Statistical analysis per policy
outcomes_df = pd.DataFrame.from_dict(outcomes3)
outcomes_df = pd.concat([experiments['policy'], outcomes_df], axis=1)
scenarios = pd.DataFrame(outcomes_df)
In [72]:
scenarios.to_csv('scenarios.csv')
In [74]:
def s_to_n(data3, direction):
mean = np.mean(data3)
std = np.std(data3)
if direction==ScalarOutcome.MAXIMIZE:
return mean/std
else:
return mean*std
In [76]:
experiments3, outcomes3 = results3
overall_scores = {}
for policy in np.unique(experiments3['policy']):
scores = {}
logical = experiments3['policy']==policy
for outcome in dike_model.outcomes:
value = outcomes3[outcome.name][logical]
sn_ratio = s_to_n(value, outcome.kind)
scores[outcome.name] = sn_ratio
overall_scores[policy] = scores
scores = pd.DataFrame.from_dict(overall_scores).T
scores
Out[76]:
In [77]:
from ema_workbench.analysis import parcoords
data4 = scores
limits4 = parcoords.get_limits(data4)
limits4.loc[0, ['Expected Annual Damage', 'Dike Investment Costs', 'RfR Investment Costs', 'Evacuation Costs',
'Expected Number of Deaths Upstream Gelderland', 'Expected Number of Deaths Downstream Gelderland', 'Expected Number of Deaths Overijssel']] = 0
paraxes = parcoords.ParallelAxes(limits4)
paraxes.plot(data4)
#paraxes.invert_axis('max_P')
plt.show()
In [78]:
def calculate_regret(data4, best):
return np.abs(best-data4)
In [79]:
# experiments, outcomes = results
overall_regret = {}
max_regret = {}
for outcome in dike_model.outcomes:
policy_column = experiments3['policy']
# create a DataFrame with all the relevent information
# i.e., policy, scenario_id, and scores
data5 = pd.DataFrame({outcome.name: outcomes3[outcome.name],
"policy":experiments3['policy'],
"scenario":experiments3['scenario']})
# reorient the data by indexing with policy and scenario id
data5 = data5.pivot(index='scenario', columns='policy')
# flatten the resulting hierarchical index resulting from
# pivoting, (might be a nicer solution possible)
data5.columns = data5.columns.get_level_values(1)
# we need to control the broadcasting.
# max returns a 1d vector across scenario id. By passing
# np.newaxis we ensure that the shape is the same as the data
# next we take the absolute value
#
# basically we take the difference of the maximum across
# the row and the actual values in the row
#
outcome_regret = (data5.max(axis=1)[:, np.newaxis] - data5).abs()
overall_regret[outcome.name] = outcome_regret
max_regret[outcome.name] = outcome_regret.max()
In [80]:
max_regret = pd.DataFrame(max_regret)
sns.heatmap(max_regret/max_regret.max(), cmap='viridis', annot=True)
plt.show()
In [81]:
colors = sns.color_palette()
data6 = max_regret
# makes it easier to identify the policy associated with each line
# in the parcoords plot
# data['policy'] = data.index.astype("float64")
limits5 = parcoords.get_limits(data6)
limits5.loc[0, ['Expected Annual Damage', 'Dike Investment Costs', 'RfR Investment Costs', 'Evacuation Costs',
'Expected Number of Deaths Upstream Gelderland', 'Expected Number of Deaths Downstream Gelderland', 'Expected Number of Deaths Overijssel']] = 0
paraxes = parcoords.ParallelAxes(limits5)
for i, (index, row) in enumerate(data6.iterrows()):
paraxes.plot(row.to_frame().T, label=str(index), color=colors[i])
paraxes.legend()
plt.show()
In [82]:
from collections import defaultdict
policy_regret = defaultdict(dict)
for key, value in overall_regret.items():
for policy in value:
policy_regret[policy][key] = value[policy]
In [84]:
# this generates a 2 by 2 axes grid, with a shared X and Y axis
# accross all plots
fig, axes = plt.subplots(ncols=4, nrows=3, figsize=(20,20),
sharey=True, sharex=True)
# to ensure easy iteration over the axes grid, we turn it
# into a list. Because there are four plots, I hard coded
# this.
axes = [axes[0,0], axes[0,1],
axes[1,0],axes[1,1],
axes[0,2],axes[0,3],
axes[1,2],axes[1,3],
axes[2,0],axes[2,1],
axes[2,3]]
# zip allows us to zip together the list of axes and the list of
# key value pairs return by items. If we iterate over this
# it returns a tuple of length 2. The first item is the ax
# the second items is the key value pair.
for ax, (policy, regret) in zip(axes, policy_regret.items()):
data7 = pd.DataFrame(regret)
# we need to scale the regret to ensure fair visual
# comparison. We can do that by divding by the maximum regret
data7 = data7/max_regret.max(axis=0)
sns.boxplot(data=data7, ax=ax)
# removes top and left hand black outline of axes
sns.despine()
# ensure we know which policy the figure is for
ax.set_title(str(policy))
plt.show()
Step 4: Scenario Discovery¶
Drop the policies and the levers¶
In [85]:
experiments_SD = experiments3.iloc[:, :19]
In [195]:
experiments_SD['Policy Number'] = experiments3['Policy Number']
In [215]:
experiments_SD['Policy Number'].unique()
Out[215]:
In [196]:
experiments_SD
Out[196]:
Step 4: Scenario Discovery¶
PRIM: Total Deaths¶
In [197]:
outcomes3['Deaths Total'] = outcomes3['Expected Number of Deaths Upstream Gelderland'] + outcomes3['Expected Number of Deaths Downstream Gelderland'] + outcomes3['Expected Number of Deaths Overijssel']
In [198]:
from ema_workbench.analysis import prim
data8 = outcomes3['Deaths Total']
y = data8 < np.percentile(data8, 25)
prim_alg = prim.Prim(experiments_SD,y, threshold=0.1)
box1 = prim_alg.find_box()
In [199]:
box1.inspect_tradeoff()
Out[199]:
In [200]:
box1.inspect(35,style='graph')
plt.show()
In [201]:
box1.show_pairs_scatter()
fig = plt.gcf()
fig.set_size_inches(20,20)
plt.show()
In [202]:
box1.select(35)
box1.show_pairs_scatter()
fig = plt.gcf()
fig.set_size_inches(20,20)
plt.show()
PRIM discovery: Annual Damage¶
In [207]:
data9 = outcomes3['Expected Annual Damage']
y = data9 < np.percentile(data9, 25)
prim_alg = prim.Prim(experiments_SD,y, threshold=0.1)
box2 = prim_alg.find_box()
In [208]:
box2.inspect_tradeoff()
Out[208]:
In [209]:
box2.inspect(35, style='graph')
plt.show()
In [210]:
box2.select(35)
box2.show_pairs_scatter()
fig = plt.gcf()
fig.set_size_inches(20,20)
plt.show()
Vizualize the results using dimensional stacking¶
In [216]:
from ema_workbench.analysis import dimensional_stacking
In [217]:
dimensional_stacking.create_pivot_plot(experiments_SD, y)
plt.show()
